function pr=prob_ReSTAT(gg_cub,bd_cub,bdtype_cub,fchs_cub,ind_cub,idx_sw_cub,mktcat_cub,pop65_cub,hhi_cub,df_cub,eerr_cub,pre,act,theta)   %For 26 parameters with loop  %###UPDATED on 02/09/2016 switching cost expressed as positive
%### UPDATED on 03/17/2017: include the ln(pr) as the expected error for
%future flows.
%### UPDATED on 05/05/2017: to measure mktdom and swc, use the quantile of
%bed size, which are mktdom*bd[1-4] and swc*bd[1-4]; include "bdtype" in
%the inputs.
%### UPDATED on 07/16/2018: replace mktdom dummy with market share.
%### UPDATED on 04/23/2020: this is the same specification but using new
%simulated future states that include obs and unobs market characteristics.
%### UPDATED on 06/22/2020: simplify by removing the for loop, hopefully
%reduce the time in fminsearch
% ### UPDATED on 07/04/2020: include HHI and totpop65, both varying by
% vendors to capture the shifting

% ev of mktdom
ev1=gg_cub.*theta(end-bdtype_cub-3);
ev1(bd_cub==0)=0;
% cost-saving per bed
vbd=theta(fchs_cub);
vbd(fchs_cub==11|fchs_cub==12)=0;
ev2=bd_cub.*vbd;
%market FEs; coefficients for market categories: theta(25, 26, 27)
%ev_mktcat=theta(mktcat_cub+24);
ev_mktcat=theta(mktcat_cub+48);  % UPDATED on 07/04/2020: the position of the parameters changes
ev_mktcat(fchs_cub==11|fchs_cub==12|mktcat_cub==0)=0;  
% sunk cost
ev3=zeros(size(ind_cub));
ev3(ind_cub>0)=theta(ind_cub(ind_cub>0)-1);
% switching cost
ev4=idx_sw_cub.*theta(end+bdtype_cub-4);
ev4(bd_cub==0)=0;
% HHI, varying by vendors
vhhi=theta(fchs_cub+12);
vhhi(fchs_cub==11|fchs_cub==12)=0;
ev_hhi=hhi_cub.*vhhi;
% totpop65, varying by vendors
vpop65=theta(fchs_cub+24);
vpop65(fchs_cub==11|fchs_cub==12)=0;
ev_pop65=pop65_cub.*vpop65;
% SUM up
%pv=sum((ev1+ev2+ev3+ev4+eerr_cub).*df_cub,2); 
eev=mean(sum((ev1+ev2+ev3+ev4+ev_mktcat+ev_hhi+ev_pop65-eerr_cub).*df_cub,2),3);
v=reshape(eev,13,[])';
expv=exp(v);
expv(pre(:,1)>1,1)=0;
idx=sub2ind(size(expv),(1:size(expv,1))',act);
pr=expv(idx)./sum(expv,2);


